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NUMERICAL  SIMULATIONS  OF  THE 
CELLULAR  STRUCTURE  OF  DETONATIONS  IN 
LIQUID  NTTROMETHANE-REGULARITY  OF 
THE  CELL  STRUCTURE 

INTRODUCTION 

Detonation  waves  were  treated  as  steady,  one-dimensional  phenomena  until  1926, 
when  Campbell  and  VVoodhead  observed  a  spinning  detonation  front  propagating  in 
a  round  tube.  This  provided  the  first  evidence  that  detonations  are  non-steady  and 
have  a  local  three-dimensional  configuration.  Since  then,  Manson,  Edwards,  Denisov 
and  Troshin,  Soloukhin,  Strehlow,  Oppenheim,  Urtiew,  and  others,  have  carried  out 
extensive  experimental  and  theoretical  studies  of  the  structure  of  gas  phase  detonations. 
Some  of  the  major  findings  of  their  work  are  summarized  in  the  texts  by  Strehlow  [l] 
and  Fickett  and  Davis  [2]. 

Numerical  simulations  of  the  detailed  structure  of  detonation  waves  began  around 
1970.  This  delay  occurred  because  such  calculations  require  computers  with  large 
memory  and  algorithms  with  low  inherent  numerical  diffusion  to  convect  compressible 
flows  where  shocks  might  develop.  The  large  amount  of  diffusion  in  early  methods 
has  the  useful  effect  of  smoothing  out  dispersion  errors  at  discontinuities.  However, 
it  also  smears  out  any  transverse  fluctuations,  which  caused  the  solutions  to  degener¬ 
ate  into  the  one-dimensional,  steady  Chapman-Jouguet  solution.  The  first  successful 
multi-dimensional  numerical  simulations  of  detonations,  carried  out  by  Taki  and  Fu- 
jiwara  [3j,  showed  that  regardless  of  the  details  of  the  triggering  disturbance,  a  fixed 
number  of  triple  points  is  formed  when  a  one-dimensional  Chapman-Jouguet  detona¬ 
tion  is  disturbed.  Oran  et  al.  [4],  starting  from  a  slightly  oblique  shock  wave  in  a 
two-dimensional  channel,  produced  a  detonation  front  with  a  single  triple  point  bounc¬ 
ing  between  the  confining  walls.  Kailasanath  et  al.  [5]  investigated  the  problem  of 
coupling  with  the  confining  walls  and  determined  a  natural  size  of  the  detonation  cell 
width  in  a  hydrogen-oxygen-argon  mixture  which  agrees  well  with  experimants. 

Manuscript  approved  April  13.  1986. 
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The  non-uniform  nature  of  detonation  waves  in  condensed  phases  was  observed  by 
Campbell  et  al.  ’6j.  Shchelkin  [7]  suggested  that  detonations  in  condensed  phases  should 
exhibit  the  same  structure  as  those  in  gases.  Urtiew  et  al.  [8]  observed  a  criss-cross 
pattern  of  diagonal  traces,  similar  to  those  inscribed  on  soot  in  gaseous  detonations, 
indented  on  the  side  wall  of  a  square  channel  in  which  a  mixture  of  nitromethane  and 
acetone  was  detonated. 

There  are  definite  similarities  in  the  stuctures  observed  in  the  condensed  and  gas 
phases.  However,  condensed- phase  detonations  differ  in  two  respects  from  gas-phase 
detonations.  First,  since  the  acoustic  impedance  of  the  walls  is  comparable  to  that  of 
the  condensed  material  they  confine,  the  wall  expands  when  a  wave  collides  with  it. 
The  wall  expansion  generates  rarefaction  waves  that  weaken  the  detonation  front  and 
make  it  oblique  to  the  wall.  In  addition,  characteristic  acoustic  and  chemical  induction 
times  are  a  few  orders  of  magnitude  less  in  condensed  phases  than  in  gases.  As  a  result, 
the  characteristic  size  of  the  detonation  ceils  is  considerably  reduced. 

In  the  present  work  we  report  the  results  of  numerical  simulations  of  the  detailed 
structure  of  detonations  in  liquid  nitromethane.  Numerical  simulations  in  condensed 
phases  are  particularly  important  since  they  can  describe  the  structure  of  the  detonation 
front  in  opaque  materials  for  which  optical  diagnostics  are  difficult  or  impossible.  The 
degree  of  numerical  resolution  required  depends  on  the  amount  of  numerical  diffusion 
inherent  in  the  finite-difference  scheme  used.  In  this  work  we  have  used  a  fourth- 
order  Flux-Corrected  Transport  algorithm  that  can  reproduce  discontinuities  with  small 
dispersion  errors  and  minimal  numerical  diffusion.  Such  an  algorithm  makes  it  possible 
to  do  two-dimensional  calculations  of  compressible  Sows  at  reasonable  costs  in  computer 
time  and  memory.  We  have  considered  the  ideal  case  of  a  planar  detonation  propagating 
in  a  narrow  slab  of  heavily-confined  nitromethane.  The  problem  modelled  therefore 
simulates  the  configuration  near  the  central  axis  of  a  wide  tube,  far  from  the  locations 
where  wall  effects  are  important.  Our  objective  is  to  show  the  formation  and  resulting 
configuration  of  multi-dimensional  detonation  structures  in  condensed  phases,  to  get 
a  preliminary  estimate  of  a  characteristic  cell  width,  and  to  investigate  the  various 
factors  controlling  the  regularity  of  the  structure. 
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PHYSICAL  MODEL 


We  have  made  the  following  assumptions: 

1.  The  detonation  is  planar,  and  thus  can  be  described  by  a  two-dimensional  model; 

2.  The  normal  component  of  velocity  at  the  walls  vanishes,  while  no  restriction  is 
imposed  on  the  tangential  velocity  component; 

3.  The  transformation  of  fuel  to  products  is  described  by  a  two-step  reaction  model 
and  occurs  in  the  bulk  of  the  condensed  fuel;  and 

4.  Within  a  small  volume  element,  the  condensed  fuel  and  the  gaseous  products  reach 
pressure  and  temperature  equilibrium  instantaneously. 

Below  we  discuss  the  implications  of  these  assumptions  on  our  calculations. 

When  a  detonation  propagates  in  a  rectangular  or  square  tube,  “transverse”  waves 
move  across  the  larger  dimension  (width)  of  the  tube  cross  section,  while  “slapping” 
waves  move  across  the  normal  direction  (height).  This  three-dimensional  configuration 
is  called  a  “rectangular”  detonation.  If  the  slapping  waves  are  suppressed,  the  front 
becomes  cylindrical.  This  two-dimensional  structure  is  called  a  “planar”  detonation  [l], 
which  is  the  subject  of  this  paper. 

Since  the  acoustic  impedance  of  a  wall  confining  a  condensed  explosive  is  compara¬ 
ble  to  that  of  the  explosive,  the  wall  expands  when  a  pressure  wave  collides  with  it.  The 
wall  expansion  generates  rarefaction  “failure”  waves  that  causes  the  detonation  front 
to  decay  near  the  wall.  Boundary  layers  also  weaken  the  structure  of  the  front  near  the 
wall.  In  our  calculations,  we  assume  that  the  normal  velocity  component  at  the  wail 
vanishes.  This  is  equivalent  to  assuming  a  perfectly  rigid  wall,  known  as  “heavy  con¬ 
finement.”  We  also  impose  no  restriction  on  the  tangential  velocity  component  at  the 
wall,  and  as  a  result  boundary  layers  cannot  form.  This  is  known  as  “free  slip.”  These 
two  assumptions  mean  that  the  walls  are  perfect  reflectors,  and  may  be  thought  of  as 
planes  of  symmetry  in  a  channel  many  times  wider.  The  resulting  detonation  structure 
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then  models  the  structure  found  in  the  inner  regions  of  a  detonating  explosive,  far  from 
the  boundaries. 

We  have  used  a  two-step  reaction  model  to  represent  the  transformation  of  fuel  to 
products.  The  first  step  is  a  chemical  induction  period,  in  which 

Fuel  — ►  Radicals,  ( step  1) 


while  in  the  second  step, 

Fuel  -r  Radicals  — ►  Products.  ( step  2) 

Since  the  fuel  we  are  considering  is  a  liquid,  the  transformation  of  fuel  to  gaseous 
products  can  occur  in  the  gas  phase,  in  the  liquid  phase,  or  at  the  interface  between 
the  two  phases.  For  example,  consider  a  small  volume  element  containing  both  liquid 
fuel  and  gaseous  products.  Heat  from  the  hot  products  could  be  transfered  to  the 
liquid  surface  causing  evaporation  of  a  thin  layer  of  fuel,  followed  by  both  reaction 
steps  occuring  in  the  gas  phase,  away  from  the  liquid  surface.  Alternatively,  heat  and 
radicals  might  diffuse  from  the  products  to  the  liquid  surface  where  the  second  step 
occurs  in  the  condensed  phase.  Finally,  both  reaction  steps  could  occur  in  the  bulk  of 
compressed  hot  liquid  in  chemical  and  thermal  isolation  from  the  gaseous  products. 

The  acoustic  relaxation  and  chemical  induction  times  in  liquid  nitromethane  are 
much  shorter  than  any  diffusion  time  scale.  We  therefore  neglect  the  effects  of  heat  and 
radicals  diffusion  from  the  hot  products  to  the  liquid  fuel  surface,  and  assume  that  both 
reaction  steps  occur  in  the  bulk  of  the  liquid  fuel.  When  a  portion  of  fuel  is  converted 
to  gaseous  products,  its  sudden  rise  in  pressure  is  communicated  through  pressure 
waves  to  the  surrounding  liquid.  .After  pressure  is  equilibrated,  the  temperature  of  the 
adiabatically  compressed  neighboring  liquid  is  generally  lower  than  that  of  the  expanded 
gases.  However,  the  liquid  reaches  this  temperature  before  any  appreciable  heating  can 
occur  by  thermal  conduction  from  the  hot  products.  In  addition,  the  chemical  induction 
time  is  short  in  the  range  of  temperatures  achieved  after  this  adiabatic  compression. 
Thus  the  liquid  fuel  forms  enough  radicals  within  its  bulk  to  start  ignition,  before 
any  significant  amount  of  radicals  diffuse  from  the  products  to  the  liquid  surface.  In 
summary,  we  neglect  the  effects  of  surface  chemical  reactions. 
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During  the  chemical  induction  period,  radicals  are  produced  by  step  1,  an  en¬ 
dothermic  process,  faster  than  they  are  consumed  by  step  2,  an  exothermic  process. 
The  concentration  of  radicals  steadily  increases  but  remains  so  small  that  the  amount 
of  fuel  and  energy  consumed,  and  the  amount  of  products  produced  during  the  induc¬ 
tion  period,  can  be  neglected.  The  system  temperature  and  the  rate  of  step  1  remains 
nearly  constant.  As  a  result,  the  rate  of  step  1  can  be  expressed  in  terms  of  measured 
induction  times,  ra(T,p),  while  reaction  step  2  can  be  assumed  to  start  only  after  this 
induction  time  is  elapsed.  This  is  particularly  useful  when  chemical  kinetic  rate  data 
are  not  available.  Induction  times  are  measured  at  constant  temperature,  T,  and  den¬ 
sity,  p.  Here,  we  use  a  superscript  o  to  indicate  that  r°  is  a  quasi-steady  function  of 
temperature,  T,  and  density,  p.  If  /  denotes  the  fraction  of  induction  time  elapsed  at 
time,  t,  then 

dt  r°(T,p)  ’  (  a) 

where  /( 0)  =  0,  while  the  rate  of  step  2  remains  zero  until  /  >  1.  If,  within  a  small 
macroscopic  volume  element,  each  of  the  radicals  can  freely  collide  with  any  of  the  fuel 
molecules,  the  rate  of  step  2  would  be  proportional  to  the  product  of  fuel  and  radicals 
concentrations.  This  is  the  case  in  gas  phase  reactions.  In  a  condensed  phase,  however, 
the  radicals  produced  are  caged  in  the  vicinity  of  their  parent  molecules,  such  that 
the  radicals  produced  at  one  location  cannot  freely  collide  with  distant  fuel  molecules. 
Therefore,  we  assume  that  the  rate  of  fuel  consumption  can  be  expressed  as 

y-{fuel  mass)  =  -(fueimass)Ae-£'/,Rr)  (2a) 

dt 

once  the  radicals  concentration  reaches  a  critical  level  defined  by  /  =  1.  Under  this 
assumption,  /  can  be  interpreted  as  the  fraction  of  critical  radical  concentration  per 
unit  mass. 
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When  the  reaction  occurs  in  a  moving  system,  d/at  in  Eqs.  (la)  and  (2a)  denotes 
a  substantial  derivative  that  follows  the  system.  Since  the  system  mass  is  always 
conserved,  if  we  denote  by  w  the  fraction  of  fuel  mass  in  a  mixture  of  fuel  and  products, 
Eq.  (2a)  can  be  written  as 

—  =  -wAe-E'RT.  (2b) 

dt 


For  fuel,  vj  =  1  yields  a  finite  rate  of  fuel  consumption  once  the  value  of  /  reaches  1. 
For  products,  w  =  0  reduces  the  rate  in  Eq.  (2b)  to  zero,  thus  ending  the  reaction 
process. 

Using  the  continuity  equation, 


dp  dpu  ,  dpv  _ 
dt  +  dx  ^  dy  ~  ’ 


Eqs.  (la)  and  (2b)  are  rewritten  as 


dpf  dpfti  dpfv 

dt  dx  ^  dy 


dpw^djw*  apw  E/ RT  (5) 

dt  dx  dy  P  ’  1  1 

respectively.  In  the  above  equations,  p  denotes  the  total  (fuel  4-  products)  density, 
while  u  and  v  denote  the  particle  velocity  components  in  the  x  and  y  directions.  Equa¬ 
tions  (4a)  and  (5)  describe  the  conservation  of  radicals  and  fuel,  respectively.  The 
left  hand  side  of  Eq.  (4a)  implies  that,  within  a  time  interval  dt,  the  radicals  coming 
through  the  surfaces  of  a  control  volume  dxdy  mix  with  the  radicals  inside,  whether 
they  were  carried  in  by  the  liquid  fuel  or  the  gaseous  products.  We  have  assumed  above, 
however,  that  the  radicals  in  the  gaseous  products  cannot  diffuse  quickly  into  the  liq¬ 
uid  and  thus  cannot  contribute  significantly  to  the  reaction  process.  Equation  (la)  is 


therefore  rewritten  as 


d(fw)  _  w 
dt  ~ 
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since  w  remains  fixed  at  w  =  1  during  the  induction  time  and  dw/dt  is  independent  of 
the  radicals  concentration  /.  Using  Eq.  (3),  Eq.  (lb)  yields 


d(pvjf)  d(p-wfu)  d(pvjfv)  _  pw 

dt  dx  dy  t°  ’ 


(46) 


where  only  the  radicals  within  the  liquid  fuel  are  mixed,  since  w  =  0  for  the  products. 
Equation  (4b)  implies  that  some  mixing  occurs  among  radicals  which  are  supposed  to 
be  caged  in  different  portions  of  liquid  fuel,  a  contradiction  to  the  model  assumptions. 
This  mixing,  however,  vanishes  as  both  dx  and  dy  tend  to  zero,  and  is  thus  part  of  the 
discretization  error- 

induction  time  measurements  in  liquid  nitromethane  were  compiled  by  Chaiken  [9], 
who  expressed  r°  in  an  Arrhenius  form,  A°exp(E°/RT).  Of  the  different  sets  of  data 
compiled,  we  selected  those  taken  at  Los  Alamos  National  Laboratory,  which  gives 
A0  =  2.3  x  10“6  pa  and  E°  =  29.1  Kcal.  The  values  of  A  and  E  in  Eq.  (2)  were  taken 
from  Mader  [lOj:  A  =  4.0  x  108  pa~1  and  E  =  53.6  Kcal. 

Finally,  we  assume  that  within  a  small  volume  element,  the  liquid  fuel  and  the 
gaseous  products  reach  pressure  and  temperature  equilibrium  instantaneously.  As  ex¬ 
plained  above,  after  pressure  equilibration  within  a  small  volume  element,  the  tempera¬ 
ture  of  the  liquid  fuel  is  generally  lower  than  that  of  the  gaseous  products.  However,  in 
order  to  avoid  the  physical  and  numerical  complications  of  solving  a  two- temperature 
problem,  we  assume  thermal  equilibrium  between  the  two  phases. 

The  set  of  conservation  equations  is  completed  by  the  momentum  equations, 


d{pu) 

a(puu) 

d(puv) 

_  _£P 

(6a) 

at 

dx 

dy 

dx  ’ 

d{pv) 

at 

a(puv) 

dx 

d(pvv) 

dy 

_ d? 

dy' 

(66) 

and  the  energy  equation, 


In  Eq.  (7),  E  =  p(e  —  {'i~  —  v2)/2),  where  e  is  the  internal  energy  per  unit  mass  and  p  is 
the  pressure.  The  conservation  relations  are  supplemented  with  the  equation  of  state. 

p  =  p(p,  e,  w) 

and 

T  =  T(p,e,w). 

The  HOM  equations  of  state,  described  by  Mader  [  10] ,  are  used  for  both  condensed 
fuel  and  gaseous  products.  The  equation  of  state  for  the  condensed  phase  is  based  on 
the  Christian  and  Walsh  technique  [llj.  This  technique  requires  measuring  the  relation 
between  the  shock  velocity  and  particle  velocity  in  the  shocked  material.  A  linear  St  to 
the  measurements  is  then  solved  simultaneously  with  the  Hugoniot  relations  to  derive 
the  conditions  behind  the  shock.  States  off  the  Hugoniot  are  related  to  those  on  it  by 
using  the  Griineisen  7  =  (dp/de) \t/p  ,  for  the  liquid. 

The  equation  of  state  of  the  gaseous  products  is  constructed  by  solving  the  equi¬ 
librium  Chapman-Jouguet  (CJ)  detonation  problem  using  the  BKW  equation  of  state 
for  the  final  products.  The  equilibrium  isentrope  through  the  CJ  point  is  evaluated 
and  the  states  off  the  isentrope  are  derived  using  Griineisen  7  for  the  products.  In 
those  regions  which  contain  a  mixture  of  fuel  and  products,  0  <  w  <  1,  the  specific 
volume  and  internal  energy  are  subdivided  between  the  two  phases  such  as  to  yield  the 
same  pressure  and  temperature.  The  details  of  the  equation  of  state  and  the  iteration 
procedure  required  to  satisfy  temperature  and  pressure  equilibrium  are  described  by 
Guirguis,  et  al.  T2i. 


S 


NUMERICAL  SOLUTION 


Equations  (3)-(7)  are  solved  by  operator  splitting  the  3uid  dynamics  and  the 
chemical  reactions  within  the  overall  computational  time  step.  For  the  fluid  dynamics, 
Eqs.  (3),  (6a),  (6b),  (7)  and  the  left  hand  sides  of  Eqs.  (4b)  and  (5)  are  solved  by 
splitting  the  d/dx  and  d/dy  operators.  Each  is  advanced  using  a  fourth-order  Flux- 
Corrected  Transport  (FCT)  algorithm  [13].  During  the  chemical  step,  the  right  hand 
side  of  Eq.  (4b)  is  advanced  using  an  explicit  Euler  method.  When  /  >  1,  indicating 
the  end  of  the  induction  time,  the  right  hand  side  of  Eq.  (5)  is  advanced  using  a 
fully- implicit  Euler  method. 

In  order  to  obtain  the  multi-dimensional  structure  of  the  detonation  wave,  the 
detailed  structure  of  the  induction  zone  has  to  be  resolved.  Otherwise,  the  solution  de¬ 
generates  into  a  one-dimensional  steady  configuration.  In  the  multi-dimensional  config¬ 
uration,  the  shock  front  is  divided  by  a  system  of  triple  points  into  alternating  stronger, 
nearly- normal,  Mach  stems,  and  weaker,  oblique,  incident  shocks.  The  reaction  front 
follows  the  Mach  stems  closely,  while  the  lower  temperatures  behind  the  incident  waves 
yield  induction  zones  extending  beyond  the  reflected  shocks,  called  here  “transverse” 
shocks.  In  order  to  obtain  the  multi-dimensional  configuration  of  the  detonation  wave, 
the  resolution  should  be  enough  to  distinguish  the  different  sizes  of  the  induction  zones 
behind  the  Mach  stem  and  incident  shock.  This  requires  keeping  several  computational 
grid  points  within  the  induction  zone  behind  the  incident  waves. 

To  estimate  the  resolution  required,  one-dimensional  calculations  were  used  as  a 
guide.  It  was  found  that  for  a  nearly-steady  one-dimensional  detonation  wave,  the 
leading  shock  wave  propagates  at  0.68  cm/^s,  while  the  temperature  behind  it  is  ap¬ 
proximately  2700  K.  The  width  of  the  induction  zone,  the  difference  between  shock 
and  particle  velocities  times  the  induction  time,  was  found  to  be  2  x  I0~4  cm.  The 
particle  velocity  is  evaluated  from  the  shock  velocity  using  the  linear  relation  between 
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the  two  velocities  described  above  in  connection  with  the  equation  of  state  of  liquid 
nitromethane. 

In  our  two-dimensional  calculations,  the  grid  along  the  y-direction  (the  width  of 
the  channel)  is  kept  fixed,  with  a  uniform  spacing  of  4  x  10~3  cm.  We  used  75,  100, 
125,  and  150  computational  cells  along  the  y-direction  giving  computational  domains 
0.03,  0.04,  0.05,  and  0.06  mm  wide,  respectively.  The  computational  grid  along  the 
x-direction  is  composed  of  200  cells,  each  initially  4  x  10“ 5  cm  wide.  When  the  shock 
front  comes  within  10  cells  from  the  right  boundary,  the  grid  expands  to  guarantee  a 
uniform  spacing  of  4  x  10“5  cm  at  all  times  around  the  shock  and  reaction  zones.  The 
result  is  that  the  shock  front  is  locked  within  cell  number  191.  The  resolution  behind 
the  induction  zone  smoothly  decreases  towards  the  left  boundary.  This  type  of  grid 
motion,  however,  limits  the  time  step  to  1/3  the  Courant  time  step,  as  explained  below. 

Denoting  the  local  grid  velocity  by  ug,  FCT  guarantees  positivity,  provided 

<  1/2,  (8) 


which  imposes  an  upper  limit  on  both  particle  and  grid  displacements  per  time  step. 
Let  us  assume 


i  l6t 

I*'i5£ol 

The  Courant  condition  for  stability  also  imposes  an  upper  limit  on  the  particle  dis¬ 
placement  in  a  single  time  step.  Let 


(]«j  +  C)g  </?<!, 


where  C  is  the  local  sound  speed  [12].  The  above  equations  have  to  be  satisfied  every¬ 
where,  and  in  particular  at  the  shock  front.  Denoting  the  front  velocity  by  U,  \ug\  >  \U\ 
at  the  shock  front.  This  is  required  for  the  right  boundary  of  the  grid  to  be  able  to 
expand  with  the  shock  and  can  be  satisfied,  if  we  choose  a  >  0,  since  jtt|  C  >  'LH 
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behind  the  shock  front  of  a  detonation  wave.  In  addition,  {it!  <  C  behind  the  shock 
front,  yielding 

,  ,  St  ^  ,3 
'“'fa  -  2- 

Using  the  above  inequalities,  Eq.  (8)  gives 


i  |  St  ,  -  St  .  .  St  ^  1  .  1  3 

|—* •'faSUfa-’-WfaS  j/9  +  «S  j«  +  <— j« 

<1. 

-  2 


Hence  a  <  1/3,  and  3  <  a  <  1/3.  In  our  calculations,  we  select  (jitj  4-  C)6t/$x  <  1/4, 
yielding  an  average  time  step  0.8  x  10-5  pts. 


RESULTS 


The  calculations  are  performed  in  a  two-dimensional  domain  representing  a  channel 
closed  at  one  end.  The  channel  is  chosen  to  be  approximately  one  detonation  cell  wide. 
The  calculations  start  with  a  one-dimensional  overdriven  detonation  which  decays  as  it 
propagates.  Before  reaching  Chapman-Jouguet  conditions,  a  finite  amount  of  energy  is 
deposited  into  a  rectangular  pocket  centered  on  the  axis  behind  the  shock-front.  When 
the  shock  wave  generated  from  the  pocket  reaches  the  detonation  front,  two  transverse 
waves  are  formed,  as  illustrated  by  the  pressure  contours  in  Fig.  1.  Figure  1  also  shows 
the  details  of  the  interaction  between  the  shock  wave  generated  by  the  pocket  and  the 
rear  wall.  After  a  few  collisions  with  the  channel  walls,  the  transverse  waves  establish 
the  front  structure. 

The  detonation  cell  width  is  estimated  from  the  calculations  by  systematically 
varying  the  width  of  the  channel.  If  the  selected  channel  is  too  narrow,  sonic  waves 
have  enough  time  to  relax  any  transverse  perturbations,  and  the  solution  degenerates 
to  a  one-dimensional  detonation.  On  the  other  hand,  if  the  channel  is  too  wide,  more 
than  one  detonation  cell  forms.  If  the  channel  is  slightly  smaller  or  larger  than  one 
cell  width,  the  sections  of  the  front  near  the  walls  intermittently  change  between  one 
and  multi-dimensional  configurations  [5].  To  start  the  trial  and  error  procedure,  the 
cell  width  was  assumed  to  be  30  times  the  induction  zone  thickness  behind  the  one¬ 
dimensional  detonation.  This  value  was  chosen  by  analogy  to  gas  phase  detonations, 
in  which  the  detonation  ceil  width  is  one  to  two  orders  of  magnitude  larger  than  the 
induction  zone  behind  a  one-dimensional  Chapman-Jouguet  detonation  [l]. 

Figure  2  shows  the  pressure  contours  behind  the  detonation  front  as  it  propagates 
in  a  0.06  mm  wide  channel.  The  pressure  contours  are  plotted  on  equally  scaled  x  and  y 
axes.  A  heavy  concentration  of  contours  occurs  at  both  shock  and  reaction  fronts.  The 
solid  traces  connecting  the  contours  at  different  times  represent  the  loci  of  the  main 
triple  points,  defined  as  those  points  on  the  front  at  which  the  shock  front  sharply 
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changes  its  curvature.  The  set  of  solid  traces  define  the  main  cellular  pattern.  The 
dashed  traces  define  a  secondary  cellular  pattern,  inscribed  by  a  set  of  secondary  triple 
points.  A  secondary  triple  point  causes  a  minor  change  in  the  shock  front  curvature 
and  is  mainly  recognized  by  the  transverse  wave  associated  with  it.  The  collection  of 
solid  and  dashed  traces  exhibit  an  irregular  cellular  pattern.  While  transverse  pressure 
waves  are  scattered  behind  those  sections  of  the  front  that  are  nearly  one-dimensional, 
they  are  localized  at  the  points  of  change  of  curvature  behind  the  most  curved  sections. 
When  two  transverse  waves  collide,  a  pair  of  triple  points  is  formed.  We  note  here 
that  the  secondary  triple  points  appears  to  originate  where  the  front  is  almost  one¬ 
dimensional. 

The  sequence  in  Fig.  2  starts  at  time  step  2500,  where  the  main  transverse  waves 
have  just  reflected  from  the  walls.  At  that  time,  the  shock  front  is  generally  symmetri¬ 
cal.  At  step  3300,  the  front  is  asymmetrical,  becoming  nearly  one-dimensional  during 
steps  3600  to  4600.  Disturbances  behind  the  front  form  new  triple  points  and  at  step 
5500,  the  channel  width  is  almost  divided  into  two  cells.  Some  of  the  structure  regu¬ 
larity  and  symmetry  is  regained  at  step  6000.  However,  during  steps  6000  to  6600,  the 
front  becomes  again  nearly  one-dimensional,  and  the  same  behavior  is  repeated.  The 
pattern  traced  by  the  triple  points  alternately  divide  the  channel  width  into  one  or  two 
cells. 

Figures  3a  and  3b  show  the  details  of  the  evolution  of  the  detonation  structure 
between  time  steps  5000-6000  and  7500-8200,  when  the  width  of  the  channel  is  divided 
into  two  cells.  In  these  two  figures  the  individual  contours  are  not  separated  to  scale, 
so  that  the  solid  and  dashed  traces  merely  connect  the  triple  points.  These  traces  are 
used  to  show  the  origin  of  the  triple  points  and  how  they  evolve.  Again,  triple  points 
are  mainly  recognized  by  changes  in  the  curvature  of  the  front  and  by  the  transverse 
waves  associated  with  them. 
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The  temperature  contours  illustrated  in  Figs.  4a  and  4b  provide  an  alternative 
approach  to  locating  the  triple  points.  The  reaction  front  is  well  defined  on  the  tem¬ 
perature  contours.  When  the  temperature  behind  the  shock  front  is  below  the  selected 
minimum  contour  level,  the  shock  front  does  not  show  on  the  contours.  This  appears 
in  the  figures  as  large  gaps  at  the  front.  The  region  confined  between  the  shock  and 
reaction  front  is  the  induction  zone.  At  a  triple  point,  there  is  a  change  in  induction 
zone  thickness  and  a  slip  line  originates.  The  stonger  the  triple  point,  the  sharper  the 
slip  line  appears  and  the  larger  the  change  in  induction  zone  thickness.  The  strength 
of  a  triple  point  is  defined  as  the  pressure  ratio  across  the  transverse  wave.  At  steps 
3200  and  3300,  in  Fig.  4a,  two  triple  points  are  colliding.  We  notice  the  sudden  change 
in  induction  zone  at  the  point  of  collision  as  well  as  the  two  slanted  slip  lines.  We  can 
also  barely  distinguish  two  transverse  shocks,  nearly  parallel  to  the  center  line  of  the 
channel.  Later  on,  at  step  3900,  we  can  distinguish  two  intersecting  slip  lines  originat¬ 
ing  from  the  two  triple  points  now  travelling  towards  the  walls.  However,  the  change 
in  induction  zone  thickness  in  this  case  appears  to  be  of  the  same  order  of  magnitude 
as  the  change  caused  by  the  disturbances  behind  the  shock  front.  As  will  be  explained 
later,  this  is  one  of  the  main  reasons  for  the  irregularity  of  the  cell  structure.  Finally, 
at  3tep  3100,  both  the  slip  lines  and  the  changes  in  induction  zone  are  barely  defined, 
although  we  can  distinguish  the  transverse  shocks  associated  with  the  two  triple  points, 
now  travelling  toward  the  center  of  the  channel. 

The  temperature  contours  in  Fig.  4b  show  a  set  of  pressure  waves  and  slip  lines 
dividing  the  field  behind  the  detonation  wave  into  many  regions.  The  two  triple  points 
near  the  center  are  in  the  process  of  colliding,  while  those  two  near  the  walls  have 
already  collided.  We  can  distinguish  two  transverse  shocks  emanating  from  the  two 
triple  points  near  the  center  and  four  slip  lines.  Following  the  sequence  from  step 
7600  to  7900,  we  see  one  of  the  transverse  shocks  interact  with  a  slip  line.  The  main 
interaction  occurs  at  step  7800.  Again,  the  change  in  induction  zone  thickness  at  the 
triple  points  is  barely  noticeable  except  when  two  triple  points  are  colliding. 
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Figures  5  and  6  show  the  pressure  contours  near  the  propagating  detonation  front 
for  two  channel  widths,  0.05  and  0.04  mm,  respectively.  Again,  the  solid  traces  represent 
the  loci  of  the  main  triple  points,  the  dashed  traces  define  a  secondary  cellular  pattern, 
and  the  collection  of  these  lines  exhibit  an  irregular  cellular  pattern. 

In  Fig.  5,  most  of  the  secondary  triple  points  are  formed  between  steps  2400  and 
2800,  when  the  front  is  nearly  one-dimensional.  The  main  cellular  pattern  is  quite 
symmetrical.  Starting  from  step  4600,  the  channel  width  is  divided  into  two  clearly 
defined  main  cells.  In  addition  to  two  triple  points  near  the  center  at  step  4600,  we 
see  a  third  main  triple  point  near  the  top  wall  where  there  is  a  change  of  the  front 
curvature.  This  also  occurs  at  step  4900,  near  the  bottom  wall. 

In  Fig.  6,  the  detonation  front  is  nearly  one-dimensional  between  steps  2100  to 
2900.  During  this  period,  the  main  triple  points  weaken  and  eventually  become  sec¬ 
ondary  triple  points,  while  three  secondary  triple  points  strengthen  and  become  main 
ones.  As  a  result,  the  main  cellular  pattern  becomes  skewed,  composed  of  one  main  cell 
reflecting  up  and  down  between  the  two  walls.  At  step  4500,  this  cell  width  measures 
approximately  0.03  mm. 

The  later  steps  in  Figs.  2,  5,  and  6  show  the  system  of  triple  points  dividing 
the  width  of  the  channel  into  more  than  one  cell,  each  approximately  0.03  mm  wide. 
However,  the  irregular  cellular  pattern  is  not  a  transient  in  the  process  of  establishing 
a  cell  width  of  0.03  mm.  This  is  verified  by  the  calculations  shown  in  Fig.  7,  which 
illustrates  the  structure  of  the  detonation  front  as  it  propagates  into  a  channel  0.03  mm 
wide.  Here  again  the  first  transverse  waves  are  produced  by  depositing  a  finite  amount 
of  energy  into  a  rectangular  pocket  centered  on  the  axis  behind  the  shock-front  of  the 
initially  one-dimensional  detonation.  Two  triple  points  are  thus  formed.  However,  after 
a  transient  period,  a  structure  is  formed  at  the  detonation  front  which  has  only  one 
main  triple  point  that  bounces  between  the  two  walls.  This  is  shown  starting  from  step 
5800  in  Fig.  7.  This  triple  point  traces  a  pattern  equivalent  to  half  a  detonation  cell. 
Since  the  walls  are  perfect  reflectors,  we  would  expect  to  observe  one  clearly  defined 
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full  cell  when  the  channel  width  is  doubled.  This  is  not  observed  however  in  a  channel 
0.06  nun  wide,  as  shown  by  Fig.  2.  We  conclude  that  the  patterns  observed  in  Fig.  2, 
5,  and  6,  are  indeed  irregular  structures,  not  just  a  transient  behavior. 
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STABILITY  OF  THE  MULTI-DIMENSIONAL  STRUCTURE 


The  formation  of  the  multi-dimensional  structure  of  a  detonation  front  can  be  un¬ 
derstood  by  considering  the  effects  of  fluctuations  on  a  homogeneously  heated  reactive 
material.  Fluctuations  in  density  or  temperature,  which  always  exist,  give  rise  to  local 
hot  regions.  In  energetic  materials,  the  chemical  induction  time  is  strongly  dependent 
on  temperature.  Once  a  hot  spot  is  formed,  it  ignites  before  the  surrounding  medium 
ignites.  Compression  waves  then  propagate  the  resulting  local  rise  in  pressure  in  all 
directions.  Because  of  the  heating  associated  with  compression,  these  waves  accelerate 
the  chemical  reaction  in  neighboring  regions.  In  turn,  these  regions  ignite  and  generate 
more  compression  waves.  Moreover,  when  the  pressure  waves  generated  from  other  hot 
spots  interact  at  some  location  with  the  first  set  of  waves,  they  create  new  hot  spot3 
there.  A  general  conclusion  is  that  truly  homogeneous  combustion  cannot  persist  in  a 
reactive  material  in  which  the  induction  time  is  a  strong  function  of  temperature. 

The  existence  of  a  one-dimensional  detonation  wave  implies  homogeneous  combus¬ 
tion  in  the  transverse  direction  behind  the  leading  shock  front.  Since  such  homogeneity 
cannot  persist  when  the  system  is  perturbed,  the  one-dimensional  configuration  is  un¬ 
stable  to  transverse  perturbations.  Erpenbeck  [14]  has  shown  that  if  the  exothermic 
kinetics  of  a  system  are  strongly  dependent  on  temperature,  one-dimensional  deto¬ 
nation  waves  are  linearly  unstable  to  transverse  perturbations.  However,  the  linear 
analysis  cannot  tell  us  what  the  final  configuration  must  be. 

When  perturbed,  the  detonation  evolves  into  a  dynamically  stable  configuration, 
a  repeatable,  organized  structure  conforming  to  the  geometry  of  the  boundaries  or 
symmetry  planes.  The  transverse  direction  is  divided  by  a  system  of  triple  points  into 
alternating  hotter  and  colder  regions.  The  regions  behind  the  nearly  normal  Mach  stems 
are  hotter.  Those  confined  between  the  incident  and  transverse  sections  of  the  front 
are  colder.  This  configuration  is  stable  because  the  induction  time  of  the  colder  regions 
is  too  long  for  fluctuations  to  start  ignition.  The  colder  regions  do  not  ignite  until  the 
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transverse  waves,  themselves  small  detonations  propagating  in  the  transverse  direction, 
consume  them.  As  the  detonation  propagates,  the  tripie  points  trace  a  cellular  pattern 
called  “cellular  structure.” 

The  multi-dimensional  configuration  described  above  is  not  always  stable,  thus  not 
always  repeatable.  As  explained  above,  one-dimensional  detonations  cannot  persist.  As 
a  result,  if  a  section  of  the  front  becomes  nearly  one-dimensional,  the  disturbances  be¬ 
hind  the  shock  front  coalesce  into  two  or  more  transverse  waves.  When  these  waves 
^  interact,  secondary  triple  points  are  formed,  leading  to  the  irregular  structures  pre¬ 

sented  above.  The  results  also  show  that  while  transverse  waves  are  scattered  behind 
those  sections  of  the  front  that  are  nearly  one-dimensional,  they  axe  localized  near  the 
triple  points  behind  the  curved  sections.  We  can  conclude  thus  that  the  sharper  the 
curvature  of  the  front  between  the  triple  points,  the  more  stable  it  is.  We  then  ask: 
what  is  the  meaning  of  increasing  the  curvature  of  the  front,  and  which  properties  of 
the  system  control  it? 

As  explained  above,  the  multi-dimensional  structure  is  stable  because  the  colder 
region  behind  an  incident  wave  does  not  ignite  until  the  reaction  front,  coming  from  the 
hotter  region  behind  the  Mach  stem,  consumes  it.  The  temperature  contours  at  steps 
3900  and  5100,  Fig.  4a,  show  that  the  change  in  induction  zone  at  the  triple  points  is  of 
the  same  order  of  magnitude  as  the  changes  caused  by  perturbations  behind  the  shock 
front.  This  observation  suggests  that  the  difference  in  induction  times,  corresponding  to 
the  difference  in  temperatures  behind  the  Mach  stems  and  incident  waves,  may  not  be 
enough  to  yield  a  stable  detonation  structure.  If  the  change  in  induction  zone  thickness 
at  the  triple  points  were  much  larger  than  any  possible  change  due  to  perturbations 
behind  the  shock,  the  multi-dimensional  structure  could  become  stable.  Below  we  show 
that  this  hypothesis  is  correct.  We  make  the  multi-dimensional  structure  more  stable 
by  increasing  the  induction  time  behind  the  incident  shocks,  the  colder  regions,  while 
decreasing  the  induction  time  behind  the  Mach  stems,  the  hotter  regions. 
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Consider  Fig.  8,  which  shows  the  Arrhenius  line  0,  the  logarithm  of  the  induc¬ 
tion  time  and  reciprocal  of  temperature.  The  horizontal  line,  included  for  reference, 
corresponds  to  the  hypothetical  case  in  which  the  induction  time  is  independent  of 
temperature.  In  this  case,  a  hot  spot  does  not  ignite  any  faster  than  the  surround¬ 
ings.  This  makes  homogeneous  combustion  and  in  turn  one-dimensional  detonations 
possible.  We  rotate  the  induction  line  to  position  1  which  is  50%  steeper.  In  order 
to  keep  the  detonation  cell  width  roughly  fixed,  we  rotate  the  induction  line  about 
a  point  on  the  line  corresponding  to  the  temperature  behind  the  shock  wave  of  the 
one-dimensional  Chapman-Jouguet  (CJ)  detonation.  Here  we  assume  that,  as  in  gas 
phase  detonations,  the  cell  width  is  proportional  to  the  induction  zone  thickness  of  the 
one-dimensional  CJ  detonation.  In  addition,  the  temperature  behind  the  shock  of  the 
one  dimensional  detonation  should  be  somewhere  in  between  the  minimum  tempera¬ 
ture  reached  behind  the  Mach  stem  and  the  maximum  temperature  reached  behind 
the  incident  shock.  These  temperatures  change  across  the  transverse  direction,  as  the 
detonation  propagates  through  one  cycle  (half  a  cell  length),  and  from  one  cycle  to  the 
next.  Rotating  the  induction  line  is  thus  equivalent  to  making  the  induction  times  be¬ 
hind  the  Mach  stems  shorter  and  behind  the  incident  waves  longer.  In  an  experiment, 
this  could  accomplished  by  mixing  the  nitromethane  with  another  liquid  such  as  to 
make  the  induction  time  a  more  sensitive  function  of  temperature. 

Figure  9  illustrates  the  pressure  contours  behind  the  detonation  front  when  the 
induction  time  is  described  by  line  1  in  Fig.  8.  Line  1  is  obtained  by  rotating  line  0 
about  2700  K,  the  temperature  behind  the  3hock  of  the  one-dimensional  detonation, 
to  a  50%  larger  slope.  Then  A0  =  1.55  x  10“ '  fis  and  E°  ~  43.65  Kcal.  All  other 
parameters  are  the  same  as  in  the  case  of  Fig.  5.  Comparing  Figs.  9  and  5  shows  that 
the  detonation  structure  has  become  much  more  regular  as  a  result  of  the  change  in 
temperature  dependence  of  the  induction  time.  The  curvature  of  the  shock  front  and 
the  two  main  triple  points  are  sharply  defined  throughout  the  whole  sequence.  We  do 
not  observe  any  secondary  triple  points.  The  structure  is  symmetrical. 
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The  details  of  the  structure  axe  shown  on  the  pressure  and  temperature  contours 
in  Figs.  10a  and  10b.  On  the  pressure  contours,  the  transverse  shocks  are  sharply 
defined,  while  a  reaction  front  appears  as  a  spotted  stripe.  The  temperature  contours 
are  heavily  concentrated  at  the  reaction  front.  On  the  temperature  contours,  the  slip 
lines  are  sharply  defined  while  the  transverse  shocks  are  barely  distinguishable.  We 
notice  the  change  in  induction  zone  thickness  near  the  triple  points. 

First  consider  the  pressure  contours  of  Fig.  10a.  The  sequence  starts  at  step  2100, 
when  the  two  triple  points  have  just  reflected  from  the  walls.  During  steps  2100  to  2900, 
the  triple  points  travel  from  the  walls  to  the  center  of  the  channel.  As  the  transverse 
shocks  interact,  from  step  2600  to  step  2900,  their  inclination  to  the  center  line  changes. 
The  curvature  of  the  shock-front  is  reversed,  and  the  Mach  stem  and  incident  shock 
switch  roles.  (The  incident  wave  is  the  section  of  the  front  that  is  squeezed  between 
the  transverse  waves).  The  reaction  front  follows  the  Mach  stem  closely,  but  is  far  from 
the  incident  shock.  On  the  temperature  contours,  two  unreacted  pockets  of  material 
appear  near  the  wall  starting  from  step  2400.  These  appear  because  of  the  formation  of 
a  second  reaction  front  behind  the  incipient  Mach  stem,  better  illustrated  in  Fig.  10b. 
Finally,  the  contours  of  step  2100  are  similar  to  those  of  step  2900,  except  for  an  upward 
or  downward  displacement  equal  to  half  the  width  of  the  channel.  The  sequence,  steps 
2100  to  2900,  describe  one  cycle  of  the  structure,  during  which  the  front  propagates 
through  half  a  cell  length. 

The  sequence  in  Fig.  10b  begins  at  a  time  just  after  the  transverse  shocks  have 
collided  and  reflected  from  each  other  at  the  center,  and  are  now  propagating  towards 
the  walls.  As  they  interact  with  the  wails  around  step  3700,  they  change  their  incli¬ 
nation,  and  by  step  3900  the  transverse  shocks  propagate  back  towards  the  center.  A 
second  reaction  front  forms  between  steps  3000  and  3300  at  the  center,  just  behind 
the  Mach  stem,  as  shown  by  the  temperature  contours.  As  the  two  reaction  fronts 
consume  the  material  between  them,  a  pocket  of  unburned  material  is  isolated  at  the 
center.  Small  pockets  of  unburned  material  should  not  affect  the  regularity  of  the 
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structure.  Large  pockets,  however,  cause  a  significant  energy- release  deficit  behind  the 
detonation,  which  eventually  disturbs  the  regularity  of  the  structure.  The  formation  of 
unreacted  pockets  has  been  observed  both  in  experiments  and  numerical  simulations 
of  gas  phase  detonations  [5, lot.  Steps  2900  to  3900  describe  a  sequence  slightly  longer 
than  one  cycle.  The  contours  of  steps  2900  and  3900  are  almost  similar,  except  for  an 
upward  or  downward  displacement  equal  to  half  the  width  of  the  channel. 

Finally,  Figs.  10a  and  10b  show  a  phase  lag  between  the  shock  front  and  the 
reaction  front  cycles.  The  transverse  shocks  always  originate  at  the  triple  points,  where 
the  shock  front  suddenly  changes  its  curvature.  However,  the  local  changes  in  the 
location  of  the  reaction  front  always  lag  behind  the  local  changes  in  the  curvature  of 
the  shock  front.  For  example,  at  step  2400,  the  reaction  front  at  the  center  is  as  close 
to  the  incident  shock  as  it  is  to  the  Mach  stem  near  the  wall.  The  reaction  front  does 
not  fall  far  behind  the  incident  shock,  as  is  expected,  until  step  2900.  At  step  2900,  the 
reaction  front  is  ahead  of  the  propagating  transverse  wave.  The  transverse  wave  does 
not  move  ahead  of  the  reaction  front  until  step  3300.  This  phase  lag  corresponds  to  a 
finite  induction  time,  since  the  reaction  occurs  because  of  the  shock  heating. 
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EFFECTS  OF  INDUCTION  TIME  PARAMETERS 


We  have  shown  in  the  previous  section  that  the  regularity  of  the  multi-dimensional 
structure  can  be  improved  by  making  the  induction  time  a  stronger  function  of  tem¬ 
perature.  This  was  accomplished  by  rotating  the  .Arrhenius  line  0  to  line  1  in  Fig.  3 
about  a  point  on  the  line  corresponding  to  the  temperature  behind  the  shock  front  of  a 
one-dimensional  detonation.  However,  rotating  the  line  also  resulted  in  the  formation 
of  large  unreacted  pockets  that  eventually  distorted  the  front  structure  again.  Below 
we  investigate  the  effects  of  rotating  the  induction  line  by  different  amounts  and  about 
different  temperatures. 

Figure  11  shows  the  results  when  the  induction  line  is  rotated  about  2700  K  to 
become  20%  steeper,  yielding  A0  =  7.33  x  10~7  /zs,  and  E°  =  34.92  Kcal.  Again,  all 
other  parameters  are  the  same  as  those  used  in  the  calculations  shown  in  Fig.  5.  The 
curvature  of  the  front  is  not  as  pronounced,  and  the  structure  is  not  as  symmetrical  as 
in  the  case  of  50%  rotation  illustrated  in  Figs.  10.  The  temperature  contours  do  not 
show  the  same  change  in  induction  zone  thickness  at  the  triple  points,  except  when 
the  two  transverse  waves  collide.  Since  the  20%  case  combines  some  of  the  features  of 
both  the  irregular  and  regular  structures,  it  was  chosen  to  test  the  convergence  of  the 
numerical  simulations. 

Figures  12a  and  12b  demonstrate  the  convergence  of  the  numerical  calculations. 
The  results  on  the  left  were  obtained  by  limiting  the  time  step  to  0.25  the  Courant  time 
step.  The  results  on  the  right  corresponds  to  the  case  when  the  time  step  was  limited 
to  0.125  the  Courant  time  step.  The  initial  conditions  of  the  second  calculations  were 
taken  as  the  resuits  of  step  2100,  from  the  drst  calculation.  The  results  are  compared 
at  nearly  equal  times.  Figure  12a  and  12b  show  the  comparison  of  the  pressure  and 
temperature  contours,  respectively. 

The  effect  of  rotating  the  induction  line  to  become  50%  steeper,  but  about  different 
temperatures,  is  illustrated  in  Figs.  13  and  14.  This  is  equivalent  to  a  combination  of 
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rotating  the  induction  line  about  2700  K,  and  a  parallel  displacement  to  shorten  or 
lengthen  the  induction  time  at  all  temperatures.  This  is  illustrated  by  lines  2  and  3 
in  Fig.  3.  All  other  parameters  are  kept  the  same  as  in  the  case  of  Fig.  2.  Figure  13 
shows  the  result  of  rotating  the  induction  line  0  about  2400  K  to  line  2,  in  which  case 
A°  =  1.11  x  10-7  ns  and  E°  =  43.65  Kcal.  the  resulting  induction  times  are  shorter 
than  those  corresponding  to  rotation  about  2700  K.  The  induction  zone  is  thinner,  thus 
reducing  the  overall  change  at  the  triple  points.  The  curvature  of  the  front  is  much 
less  than  those  curvatures  observed  in  Figs.  10,  and  the  structure  not  as  symmetrical. 
In  addition,  the  cell  appears  to  be  smaller  than  the  channel  width. 

The  results  of  rotating  the  induction  line  about  2800  K  to  become  50%  steeper 
are  shown  in  Figs.  14.  In  this  case  A°  —  1.71  x  10~7  n&  and  &0  =  43.65  Kcal, 
denoted  as  line  3  on  Fig.  8.  The  resulting  induction  times  are  slightly  longer  than 
those  corresponding  to  rotation  about  2700  K.  The  induction  zone  is  thicker,  thus 
giving  a  larger  change  in  induction  zone  thickness  at  the  triple  points.  Such  a  large 
induction  zone  increases  the  tendency  to  form  unburned  pockets,  as  illustrated  by  the 
sequence  of  temperature  contours  in  Fig.  14a.  At  step  1500,  following  collision  of  the 
triple  points  at  the  center,  a  second  reaction  front  is  formed  behind  the  new  Mach 
stem.  As  a  consequence,  a  pocket  of  unreacted  material  is  isolated  at  step  1700.  The 
pressure  contours  in  Fig.  14b  exhibit  as  sharp  a  curvature  of  the  front,  as  in  the  case 
of  rotation  about  2700  K  shown  in  Figs.  10.  However,  because  the  unburned  pocket 
holds  energy  and  releases  it  later,  the  transverse  shocks  are  slightly  dispersed  and  the 
structure  is  not  completely  symmetrical.  In  addition,  two  new  large  pockets  form  near 
the  wall.  This  is  shown  in  the  temperature  contours  starting  from  step  3100. 

As  a  compromise  between  a  large  change  of  induction  zone  at  the  triple  point 
and  the  tendency  to  form  unburned  pockets,  Figs.  15a  and  15b  show  the  results  of 
rotating  the  induction  line  to  become  only  30%  steeper,  about  a  point  on  the  line 
corresponding  to  2700  K.  In  this  case,  A0  =  4.57  x  10-7  /zs  and  E°  =  37.83  Kcal. 
All  other  parameters  in  the  calculation  are  the  same  as  those  used  in  the  calculation 


shown  in  Fig.  5.  Figure  15a  shows  in  detail  the  process  of  collision  of  the  triple  points 
with  the  walls,  while  Fig.  15b  shows  the  details  of  collision  at  the  center  of  the  channel. 
Since  the  wails  are  perfect  reflectors,  these  two  processes  should  be  similar  provided 
the  structure  is  stable  in  the  sense  that  it  is  repeatable  from  cycle  to  cycle.  The  change 
In  induction  zone  at  the  triple  points  is  not  as  large  as  in  the  case  of  the  50 %  steeper 
induction  line  shown  in  Figs.  10.  However,  no  large  unburned  pockets  are  formed. 
Thus  although  the  curvature  of  the  front  is  not  as  pronounced,  the  structure  remains 
symmetrical  and  the  transverse  shocks  are  sharply  defined  at  all  times.  We  can  clearly 
distinguish  the  slip  lines  originating  at  the  triple  points  on  the  temperature  contours. 
Finally,  the  repetition  from  cycle  to  cycle  is  improved.  Contours  which  are  1000  steps 
apart  are  very  similar,  differing  only  by  an  upward  or  downward  displacement  equal  to 
half  the  channel  width. 


DISCUSSION  AND  CONCLUSIONS 


We  have  used  time- dependent  two-dimensional  numerical  simulations  to  study  the 
detailed  structure  of  planar  detonation  waves  in  liquid  nitromethane.  To  calculate 
the  multi-dimensional  structure  of  the  front,  a  fourth-order  Flux- Corrected  Transport 
(FCT)  algorithm  was  used  for  convection.  In  addition,  it  was  necessary  to  resolve  the 
reaction  zone  with  a  number  of  grid  points.  Convergence  was  achieved  using  4x  10~5  cm 
computational  cells  and  an  average  time  step  of  0.3  x  10" 5  /xs. 

Since  we  have  neglected  the  boundary  layer  effects  and  have  assumed  rigid  wails, 
the  simulation  models  the  inner  region  of  a  large  slab  of  material.  Failure  waves, 
originating  when  the  wall  expands,  and  boundary  layers  distort  the  structure  of  the 
detonation  front  and  those  cells  near  the  wall.  Since  a  typical  detonation  cell  in  liquid 
nitromethane  is  about  0.05  mm,  any  practical  system  is  at  least  20  to  a  100  ceils  wide. 
The  inner  region  of  such  a  system  is  expected  to  contain  a  large  number  of  cells  for  which 
the  wall  effects  are  insignificant  and  therefore  can  be  described  by  heavy  confinement. 

The  simulations  show  that  when  the  one-dimensional  detonation  front  is  perturbed, 
it  evolves  into  a  system  of  colliding  triple  points  that  trace  a  cellular  pattern.  In 
pure  liquid  nitromethane,  the  resulting  cellular  pattern  is  irregular.  The  width  of  a 
detonation  cell  is  approximately  0.05  mm.  It  is  important  to  note  however  that  the 
accuracy  of  this  estimate  is  contingent  on  how  accurately  the  equation  of  state  and 
chemistry  are  represented. 

The  equation  of  state  has  been  extensively  used  in  the  past.  In  one-dimensional 
systems  it  gave  results  reasonably  close  to  experimental  observations  [10J.  The  two-step 
reaction  chemistry  model  that  we  used  is  based  on  parameters  obtained  from  direct 
measurements  of  the  induction  times  in  liquid  nitromethane.  Since  the  induction  time 
plays  the  major  role  in  forming  the  cellular  structure,  a  simple  model  that  uses  directly 
measured  induction  times  should  be  superior  to  any  detailed  chemical  kinetics  model 
in  which  the  reactions  rates  are  at  best  a  very  rough  guess. 
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Experimental  observations  of  the  cellular  structure  in  liquid  nitromethane  were 
carried  out  using  mixtures  of  nitromethane  and  acetone  [1,2, Si,  thus  the  observed  det¬ 
onation  cells  size  cannot  be  directly  compared  to  the  model  predictions.  However, 
the  cellular  structures  observed  in  mixtures  of  nitromethane  and  acetone  are  irregular, 
which  is  consistent  with  the  model  predictions. 

The  calculations  of  the  irregular  structures  showed  a  correlation  between  the  cur¬ 
vature  of  the  front  and  generation  of  disturbances  that  caused  the  irregularity.  As  a 
detonation  propagates,  its  front  often  becomes  nearly  one-dimensional.  The  simula¬ 
tions  then  show  weak  transverse  pressure  waves  scattered  behind  the  flat  sections  of 
the  front.  Coalescence  of  two  or  more  of  these  waves  then  generates  the  secondary 
triple  points  that  cause  the  non-uniform  pattern.  In  contrast  to  this,  the  transverse 
pressure  waves  are  relatively  localized  at  the  points  of  change  of  curvature  behind  the 
highly  curved  sections  of  the  front.  We  conclude  that  highly  curved  detonation  fronts 
lead  to  more  regular  patterns. 

The  calculations  of  the  irregular  structure  also  show  that  the  change  in  the  thick¬ 
ness  of  the  induction  zone  at  the  triple  points  is  of  the  same  order  of  magnitude  as  the 
changes  caused  by  the  tranverse  waves  scattered  behind  the  front.  This  led  us  to  the 
hypothesis  that  the  structures  would  be  more  regular  if  the  change  in  induction  zone 
thickness  at  the  triple  point  were  made  larger  than  any  changes  caused  by  perturba¬ 
tions  behind  the  front.  This  was  verified  by  changing  the  induction  time  so  that  it 
became  a  stonger  function  of  temperature,  and  repeating  the  calculations.  A  stronger 
function  means  a  larger  change  in  induction  time  for  a  given  change  in  temperature. 

The  degree  of  regularity  of  the  detonation  cell  structure  is  sensitive  to  the  tem¬ 
perature  dependence  of  the  induction  time.  This  was  shown  by  rotating  the  .Arrhenius 
line  (the  relation  between  the  logarithm  of  the  induction  time  and  reciprocal  of  tem¬ 
perature,  Fig.  3)  to  various  slopes  and  about  various  points  on  the  line.  The  most 
regular  structure  was  obtained  when  the  Arrhenius  line  was  made  30%  steeper  and  was 
rotated  about  a  point  corresponding  to  the  temperature  behind  the  shock  front  of  the 
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one-dimensional  detonation.  Smaller  slopes  did  not  produce  stable  multi-dimesional 
structures.  Larger  slopes  resulted  in  the  formation  of  large  pockets  of  unreacted  ma¬ 
terial,  which  eventually  distorted  the  symmetry  and  uniformity  of  the  structure  from 
one  detonation  cell  to  the  next.  In  general,  the  calculations  in  which  the  induction  line 
was  rotated  about  a  lower  temperature  showed  an  irregular  structure.  Those  in  which 
the  induction  line  was  rotated  about  a  higher  temperature  showed  a  greater  tendency 
to  form  pockets  of  unreacted  material. 

When  the  calculations  showed  a  regular  cellular  structure,  we  did  not  see  the  weak 
triple  points  that  formed  the  cellular  substructure  in  the  irregular  cases.  Only  sharply 
defined  triple  points  traced  the  cellular  pattern.  The  curvature  of  the  front  was  also 
sharply  defined  at  all  times.  The  change  in  induction  zone  thickness  at  the  triple  points 
was  much  larger  than  any  changes  caused  by  disturbances  behind  the  front.  However, 
even  when  the  structure  was  regular  it  is  not  nearly  as  repeatable  as  those  patterns 
observed  and  calculated  for  H2-O2  mixtures  highly  diluted  in  argon  [5|. 

Finally,  it  is  important  to  note  that  in  the  liquid  nitromethane  system  we  have 
studied,  the  chemical  induction  time  is  much  longer  than  the  energy  release  time. 
Specifically,  the  induction  time  is  one  order  of  magnitude,  larger  than  the  energy  re¬ 
lease  time  at  2700  K,  the  temperature  behind  the  shock  front  of  the  one-dimensional 
detonation.  Future  work  should  include  parametric  studies  of  other  effects  such  as  en¬ 
ergy  release  process  and  the  equation  of  state  on  the  regularity  of  the  cellular  structure. 
It  is  important  to  investigate  the  effects  of  having  a  channel  width  which  is  two  detona¬ 
tion  cells  or  wider.  .Another  important  approach  is  to  calculate  the  cellular  structure  in 
H2-O2  mixtures  without  any  dilution.  Kailasanath,  et  al.  [5]  obtained  highly  regular 
structures  for  mixtures  of  H2-O2  diluted  in  argon.  From  what  we  have  now  learned, 
we  expect  to  see  irregular  structures  in  H2-O2  mixtures  if  no  argon  is  present. 
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5  —  Pressure  contours  of  a  detonation  wave  propagating  in  a  channel  0.05  mm  wide.  Solid 
traces  are  loci  of  main  triple  points.  Dashed  traces  are  loci  of  secondary  triple  points. 
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10  —  Pressure  and  temperature  contours  behind  the  detonation  wave  in  Fig.  9.  (a)  Collision  of  the  two 
triple  points  at  the  center  of  channel;  (b)  Collision  of  the  two  triple  points  with  the  walls. 
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Fig.  10  (Continued)  —  Pressure  and  temperature  contours  behind  the  detonation  wave  in  Fig.  9.  <a) 
Collision  of  the  two  triple  points  at  the  center  of  channel;  (b)  Collision  of  the  two  triple  points  with  the 
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Fig.  11  —  Pressure  and  temperature  contours  behind  a  detonation  wave  propagating  in  a  channel  0.05 
mm  wide.  Induction  time  is  obtained  by  rotating  line  0  in  Fig.  8  about  2700  K  to  a  20%  larger  slope: 
*  -'.83x10'  exp(34920/«f)  Ms. 
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Fig.  12  (Continued)  —  Numerical  convergence  test  for  the  same  conditions  of  Fig.  11,  except  that 
calculation  on  the  right  start  at  step  2100  using  a  time  step  half  of  that  used  in  calculating  the  results  on 
the  left.  Comparisons  are  at  nearly  equal  times,  given  at  the  bottom  right  of  each  frame,  (a)  Pressure 
contours:  (b)  Temperature  contours. 
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(a) 

Fig.  14  —  Pressure  and  temperature  contours  behind  a  detonation  wave  propagating  in  a  channel  0.06 
mm  wide.  Induction  time  corresponds  to  line  3  in  Fig.  8,  obtained  by  rotating  tine  0  about  2800  K  to  a 
50%  larger  slope;  1.71  x  10“'  exp(43650//?D  /is.  (a)  Formation  of  a  pocket  of  unreacted 
material  at  the  center  of  channel;  (b)  Distortion  of  structure  due  to  formation  of  large  pockets  of 
unreacted  material. 
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Fig.  14  (Continued)  —  Pressure  and  temperature  contours  behind  a  detonation  wave  propagating  in  a 
channel  0.06  mm  wide.  Induction  time  corresponds  to  line  3  in  Fig.  8.  obtained  by  rotating  line  0 
about  2800  K  to  a  50%  larger  slope;  r®—  1.71  x  10~7  exp(43650//?D  /j.s.  (a)  Formation  of  a  pocket 
of  unreacted  material  at  the  center  of  channel;  (b)  Distortion  of  structure  due  to  formation  of  large 
pockets  of  unreacted  material. 
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Fig.  15  —  Pressure  and  temperature  contours  behind  a  detonation  wave  propagating  in  a  channel  0.05 
mm  wtde.  Induction  time  is  obtained  by  rotating  line  0  in  Fig.  8  about  2700  K  to  a  50%  larger  slope: 

4.57  x  10“'  exp(37830//?T)  j*s.  (a)  Collision  of  the  two  triple  points  with  the  wails:  (b)  Collision 
of  the  two  triple  points  at  the  center  of  channel. 


